%T1ScanExperiment( loadStr, saveStr, method,sliceNum)

function T1ScanExperiment(...
    loadStr, saveStr, method,sliceNum)

%Global setting
savefitdata = 1;
% Load the data 
load(loadStr);
data = IMG;

%%
% f=figure, imshow(abs(I.M(:,:,fix(size(I.M,3)/2))),[0 1000]);
% [x y]=ginput(2);
% x1 =fix(min(x));
% y1 =fix(min(y));
% x2 =fix(max(x));
% y2=fix(max(y));
% data =abs(I.M(y1:y2,x1:x2,:));
% f=figure, imshow(data(:,:,end));
% data =abs(I.M(:,:,:));

data = double(data);
[dim1, dim2,dim3] = size( data );
data = reshape( data, dim1,dim2,sliceNum,dim3/sliceNum);
% construct the extra struct
extra = struct('tVec',[200:1000:4200 300:1000:2300 ],'T1Vec',1:3000);
switch method
    case{'SASHA'}%1
        extra.tVec = Info.TIs;
        nlsS = getNLSStruct(extra,1);
        nbrOfFitParams = 6; % Number of output arguments for the fit
        fitStr = ['[T1Est bMagEst aMagEst res sigma r2] = SASHA(data(jj,:),nlsS,3);']    % 2 or 3 for fitting num
        storeStr = 'T1(jj,:) = [T1Est bMagEst aMagEst res sigma r2];'
        clearStr = 'clear jj T1Est bMagEst aMagEst res sigma r2'
    case{'MOLLI'}%2
        nbrOfFitParams = 4; % Number of output arguments for the fit
        extra.tVec = Info.Tinv;
        nlsS = getNLSStruct(extra,1);
        fitStr = ['[T1Est bMagEst aMagEst res] = MOLLI(abs(data(jj,:)),nlsS);']
        storeStr = 'T1(jj,:) = [T1Est bMagEst aMagEst res];'
        clearStr = 'clear jj T1Est bMagEst aMagEst res'
    case{'SATCOT1T2Fitting'} %3
           
        nlsS.tdVec =  Info.Tsat;
        nlsS.teVec =  Info.Tes+eps;
        nlsS.t2Wflg = [0, 0, 0, 1, 0, 0, 1, 0, 0, 1];
        nlsS.t1Wflg = [1, 1, 1, 0, 1, 1, 0, 1,1, 0 ];    
        nlsS.T1Init = 1500 ;
        nlsS.T2Init = 40 ;
        nlsS.x0 = [ 2000  1500 40 ];
        nbrOfFitParams = 5; % Number of output arguments for the fit
        fitStr = '[T1Est, T2Est, aEst, res,sd] = SATCOT1T2Fitting(data(jj,:),nlsS,1);' ;
        storeStr = 'T1(jj,:) = [T1Est, T2Est, aEst, res,sd];';
        clearStr = 'clear jj T1Est T2Est aEst res sd';
   
    case{'COSATIRT1Fitting'}   %4
          try
            nlsS.FailedMocoBit = Info.FailedMocoBit;
          catch
            nlsS.FailedMocoBit = zeros(1,10);
          end
        nlsS.tINVVec=Info.Tinv(nlsS.FailedMocoBit==0);
        nlsS.tSATVec =Info.Tsat(nlsS.FailedMocoBit==0);   
        dataPolar =   [ 0 1 1 1 1 0 0 0 0 0  ] ;
        nlsS.dataPolar = dataPolar(nlsS.FailedMocoBit==0);
        nlsS.T1Init = 1400 ;
        nbrOfFitParams = 5;
        fitStr =[ 'nlsS.x0 = [ max(data(jj,:))   nlsS.T1Init ]; '...
            '[T1Est, aEst, bEst, res,sd] = COSATIRT1Fitting(data(jj,:),nlsS,1);' ];
        storeStr = 'T1(jj,:) = [T1Est, aEst, bEst, res,sd];';
        clearStr = 'clear jj T1Est aEst bEst res sd';           
end

%%
dims = size(data);
nbrow = size(data,1);
nbcol = size(data,2);

clear tmpData;
if numel(dims) > 3
    nbslice = dims(3); % Check number of slices
else
    nbslice = 1;
    tmpData(:,:,1,:) = data; % Make data a 4-D array regardless of number of slices
    data = tmpData;
    clear tmpData;
end

dataOriginal = data();
% Mask the background (points dimmer than maskFactor*the brightest point)
% Done on the last image, where we expect the magnetization to have almost
% fully recovered.
% Increase mF is the mask does not cover all the background. Decrease it if
% it covers part of the object.
mF = 0.01;
maskFactor = mF;
mask = zeros(nbrow, nbcol, nbslice);

[u,v] = max(extra.tVec);
v = 1;
for kk = 1:nbslice
    maskTmp = mask(:,:,kk);
    maskTmp = medfilt2(maskTmp); % remove salt and pepper noise
    maskThreshold = maskFactor*max(max(abs(data(:,:,kk,v))));
    maskTmp(find(abs(data(:,:,kk,v))> maskThreshold)) = 1;
    mask(:,:,kk) = maskTmp;
    clear maskTmp
end

maskInds = find(mask);

nVoxAll = length(maskInds);
% How many voxels to process before printing out status data
numVoxelsPerUpdate = min(floor(nVoxAll/10),1000);

T1 = zeros(nVoxAll, nbrOfFitParams);
% Number of status reports
nSteps = ceil(nVoxAll/numVoxelsPerUpdate);

for ii = 1:size(data,4)
    tmpVol = data(:,:,:,ii);
    tmpData(:, ii) = tmpVol(maskInds)';
end
clear tmpVol;
data = tmpData;
clear tmpData;

startTime = cputime;
fprintf('Processing %d voxels.\n', nVoxAll);

h = waitbar(0, sprintf('Processing %d voxels', nVoxAll));

for ii = 1:nSteps
    curInd = (ii-1)*numVoxelsPerUpdate+1;
    endInd = min(curInd+numVoxelsPerUpdate,nVoxAll);
    for jj = curInd:endInd
        % Do the fit
        eval(fitStr)
        % Store the data
        try
            eval(storeStr);
        catch
            T1Est, aEst, res,sd,ir;
        end
    end
    waitbar(ii/nSteps, double(h), sprintf('Processing %d voxels, %g percent done...\n',nVoxAll,round(endInd/nVoxAll*100)));
end
eval(clearStr)
close(h);
timeTaken = round(cputime - startTime);
fprintf('Processed %d voxels in %g seconds.\n',nVoxAll, timeTaken);

dims = [size(mask) 4];
im = zeros(size(mask));

for ii = 1:nbrOfFitParams
    im(maskInds) = T1(:,ii);
    T1_temp(:,:,:,ii) = im;
end
clear T1;
T1 = T1_temp;
Mag =data;
% Going back from a numVoxels x 4 array to nbrow x nbcol x nbslice
if (savefitdata)
    save(saveStr,'T1','mask','nlsS','Mag') ;
end
